1 Inleiding

De volksbank is het moederbedrijf van ASN, SNS, BLG en Regiobank. Als de vierde bank van Nederland, financiert de volksbank hypotheken, beheert ze spaargeld en biedt de bank klanten een betaalrekening. Binnen de afdeling Financial Markets van de Volksbank, hebben de macro economen van de Volksbank als taak de toestand in de economie te monitoren en toe te lichten aan directie en beleidsmakers. Een belangrijk onderdeel daarvan heeft betrekking met de ontwikkeling van de hypotheekmarkt in Nederland. De huizenmarkt is een belangrijke factor in de beleidsontwikkeling van de bank, met betrekking tot de hypotheekproductie.

De huizenprijzen en aantallen verkochte woningen zijn de belangrijkste indicatoren voor de ontwikkeling van de hypotheekmarkt.

Voor deze analyse is de PPDAC cyclus ( Problem, Plan, Data, Analysis, and Conclusions) gevolgd. Een methode voor het uitvoeren van statististch onderzoek [OldfordR.J.MacKay2000]. De opbouw van dit verslag volgt deze cyclus. In de volgende paragraaf wordt eerst de probleemstelling nader toegelicht. Vervolgens is in de planningsfase een literatuurstudie uitgevoerd om de belangrijkste factoren te bepalen die nodig zijn om de vraagstelling te beantwoorden. De data benodigd voor dit onderzoek is daarna verzameld en vervolgens geanalyseerd. Op basis daarvan worden tot slot een aantal conclusies getrokken. Deze cyclus is een aantal keer doorlopen - op basis van de uitkomsten. Dit komt terug in de probleemstelling en de uitwerking van deze casus.

2 Probleemstelling

Uit het overleg met de macro economen is de vraag naar voren gekomen of een beter voorspellend model kan worden ontwikkeld voor de huizenmarkt in Nederland.

Is het mogelijk om modellen te ontwikkelen die de huizenprijs op middellange termijn voorspelt en het aantal transacties op de huizenmarkt op middellange termijn voorspelt?"

Vanuit de Volksbank is in het verleden al eerder statistisch onderzoek gedaan naar deze vebanden. Data werd hierbij verzameld via Datascope (Refinitiv) en geanalyseerd in Excel.

Voor de huizenprijzen is eerder gekeken naar een verband tussen huizenprijzen en de werkgelegenheid, lonen, GDP en rente. Hieruit is een eenvoudig model ontwikkeld waarbij met name het verband met GDP wordt gebruikt om de ontwikkeling van de huizenprijzen te voorspellen. Voor het aantal transacties op de lange termijn wordt een verband verondersteld met het aantal huishoudens. op kortere termijn spelen economische groei, rentestanden en evt. regelgeving ook een rol.

In deze exploratie wil ik kijken naar mogelijke factoren die effect hebben op de marktprijs. Het is in deze fase nog niet de ambitie om ook daadwerkelijk een alternatief model te ontwikkelen voor het voorspellen van de marktprijzen en aantal transacties.

Welke factoren spelen een belangrijke rol in het bepalen van de huizenprijzen en transactie volumes op een termijn van 5 tot 10 jaar?

Aanvullend zijn er een aantal vragen geformuleerd waar mogelijk door middel van exploratie iets meer over te weten kunnen komen:

  1. Nederland heeft een van de hoogste persoonlijke schulden (hypotheeklast) in de wereld. Wat zijn de gevolgen voor de hypotheekmarkt van het gewijzigde beleid? Heeft dit effect op de huizenprijzen?

  2. Hoe ontwikkeld de krapte op de huizenmarkt zich, op termijn? Momenteel zijn er ambiteuze plannen om ergens tussen de 75.000 en 100.000 woningen per jaar te gaan bouwen. Het aanbod aan huizen wordt dus groter - terwijl op termijn verwacht wordt dat de bevolking weer gaat krimpen. Kunnen we hier een voorspelling van maken?

  3. Hoe zit het met de bevolkingsgroei? Welk deel van de aanwas in Nederland komt door migratie? Hoe zit dat regionaal? Is er nog sprake van verstedelijking na corona? Stroomt het platteland leeg - en de randstad vol? Welk effect heeft dit op de huizenprijzen?

  4. Zien we demografisch verdeeld verschillen waar mensen willen wonen? Verhuizen gezinnen met kinderen naar het platteland, en ouderen weer naar de stad? Het aantal huishoudens dat uit 1 persoon bestaat is erg gegroeid in de loop der tijd - en zal (denk ik) nog wel doorstijgen met de vergrijzing. Betekent dit dat er meer mensen naar de stad willen verhuizen?

  5. Hoe afhankelijk is de huizenmarkt van de stand van de economie? Gaan huizenprijzen alleen omhoog tijdens de groeifase van de economie?

De eerste deelvraag van dit onderzoek gaat in op de vraag in hoeverre de huizenprijzen wordt bepaald door vraag en aanbod.

Daarna wordt verder ingegaan op de factoren die op langere termijn de huizenprijzen beinvloeden en wordt gekeken of hier voorspellende waarde in zit.

Vervolgens wordt gekeken naar de transactievolumes. Transactievolumes hebben waarschijnlijk een direct verband met het aantal huishoudens. Hierbij is gekeken of er seizoensinvloeden of trendmatigheden in de transactie volumes te ontdekken zijn.

Ten slot wordt kort ingegaan op de vraag of er regionaal verschillen zijn in de ontwikkeling van de huizenprijzen.

3 Plan

Om in ieder geval een begin te maken met het beantwoorden van deze vragen ben ik in gesprek gegaan met de specialisten bij de Volksbank. Aanvullend heb ik een korte literatuurstudie gedaan naar de factoren die op lange termijn de prijsstelling zouden beinvloeden.

In eerste instantie richt het onderzoek zich op de vraag en aanbod op de woningmarkt in Nederland. Het CBS publiceert vrij veel over de prijsonwikkeling van bestaande woningen, nieuwe woningen en de stand van de economie. De Socrates modellen sluiten goed aan bij deze vraagstelling (zie Kenneth Gopal, Gerard van Leeuwen, David Omtzigt, 2020).

Op langere termijn is de woningmarkt een markt van vraag en aanbod. Vraag en aanbod worden bepaald door:

Vraag woningen:

  • Aantal huishoudens
  • Demografie (omvang en samenstelling huishoudens)
  • (Arbeids)migratie

 

Aanbod woningen:

  • Woningvoorraad
  • Nieuwbouw / Sloop
  • Bouwkosten
  • Beschikbaarheid bouwlocaties
  • Ruimtelijke ordeningbeleid

Daarna is verder ingegaan op de factoren die van invloed zijn bij de prijsonwikkeling op de woningmarkt. Hierbij kan een onderscheidt gemaakt worden in factoren die direct de prijsontwikkeling beinvloeden - en factoren die meer aanduiden of er sprake is van oververhitting (Groot, Vogt, Van der Wiel, & Van Dijk, 2018).

Indicatoren

  • Gem. huizenprijzen (prijsindex)
  • Besteedbaar Inkomen
  • Particulier vermogen
  • Werkloosheid
  • Rentestanden
  • Financiering hypotheek of eigen geld
  • Fiscale behandeling, subsidies
  • Toegang tot krediet
  • Voorkeuren van kopers
  • Percentage huur / koop woningen
  • Consumentenvertrouwen

Oververhitting

  • Gem. huurprijzen (berekenen verhouding koop/huur)
  • Aantal woningen dat te koop staat
  • Gem. duur dat een woning te koop staat
  • Verschil tussen aanbodprijs en transactieprijs
  • Aanbod huizen

Het aantal transacties wordt met name bepaald door het aantal huishoudens. Daarnaast kan door krapte op de woningmarkt een schaarste ontstaan waardoor de woningmarkt ‘op slot’ geraakt. Tijdelijk zou prijsinelastisiteit er voor kunnen zorgen dat consumenten geen nieuwe woning betrekken omdat zij de huizen te duur vinden. Het consumentenvertrouwen speelt hierbij een belangrijke rol.

Indicatoren aantal transacties

4 Data

De gewenste informatie is opgevraagd bij de statistici (via Reuters Datastream) en later aangevuld met open data van het CBS (Statline).

Voor het inlezen van de (spreadsheet) Reutersdata zijn een aantal functies ontwikkeld. Om de data van het CBS gemakkelijker in te kunnen lezen is een routine ontwikkeld die aan de hand van metadata de gevraagde data inleest. Routines voor het lezen en opschonen van de data, zijn opgeslagen in seperate R files.

De data die via die routines is gegenereerd wordt tussentijds opgeslagen in Rdata formaat - en kan voor de analyse hier direct worden ingelezen. Op deze wijze is een scheiding aangebracht tussen de data-acquisitie en de analyse van de resultaten.

Ter beperking van de omvang van dit onderzoek is een selecties gemaakt op basis van de beschikbare data.

Voor de indicatoren over de huizenprijzen zal ik mij alleen beperken tot het besteedbaar inkomen, de rente standen, particulier vermogen en werkloosheid en het consumenenvertrouwen.

Gegevens over de verhouding tussen huur en koopprijzen, het aantal te koop staande huizen, de gem. tijd dat een woning te koop staat en de prijsverschillen tussen aanbodprijs en transactieprijs zeggen iets over de mate van oververhitting van de markt (zie Groot et al. (2018)). Ik beschik echter niet over voldoende data om hier een tijdsanalyse over uit te voeren.

Voor de indicatoren over het aantal transacties beperk ik mij tot een tijdreeksanalyse van de transacties zelf. Langere termijn afhankelijkheden van de economie, het consumentenvertrouwen en het aantal huishoudens zullen er zeker zijn maar zijn vanwege de beperkte tijd voor dit onderzoek niet verder meegenomen.

library (tidyverse)          # standard tidyverse packages
library (here)               # reletive paths for files
library (lubridate)          # Work with dates
library (reshape2)           # for easy melting wide format data
library (knitr)              # for outputting data tables
library (kableExtra)         # Nice layout options for data tables
library (viridis)            # Nicer colors for plots

library (forcats)            # for handling categorical data

library (forecast)           # Time series forecasting
library (timetk)             # Time series forecasting
library (fpp2)               # Time series forecasting
library (zoo)                # handling time series
library (urca)               # for time serries analyses

library (ggpubr)             # Mix multiple graphs
library (cowplot)            # for alignin multiple plots 

library(sjPlot)              # For printing model output
library(sf)                  # spacial data for geographicals maps


library(rsample)             # data splitting for model fitting
library(earth)               # fit MARS models
library(caret)               # automating the tuning process
# install.packages("vip")
library(vip)                 # variable importance
# install.packages("pdp")
library(pdp)                 # variable relationships
library(timetk)
library(corrplot)            # Correlation plot
# Loading the datasets

dfs<-list.files(here("data","tidy"), pattern = "*.Rdata")
for(d in dfs) {
  load(file= here("data","tidy",d ))
}

# Set default theme for all plots

theme_set(theme_minimal())

5 Analyse

5.1 Prijsontwikkeling woningmarkt

Prijzen van koopwoningen zijn de laatste jaren erg sterk toegenomen. We zien kleine verschillen per type woning. Waarbij met name appartementen de laatste tijd releatief duurder zijn geworden.

prijsindex_type_woning %>%
  filter(type_woning != "Totaal woningen") %>%
  filter(type_woning != "EGW") %>%
  ggplot(aes(x=date, y=prijsindex_type_woning)) +
  geom_line(aes(color = type_woning)) +
  facet_wrap(~ type_woning) +
  labs(
    title = "Prijsindex Bestaande Koopwoningen (PBK)",
    y = "",
    x = '',
    color = 'Type woning'
    ) +
  theme(legend.position = "none")

5.1.1 Vraag naar woningen

Bij de factoren die de vraag naar woningen bepalen kijken we naar de ontwikkeling van het aantal huishoudens en de samenstelling van de huishoudens. Nederlanders worden steeds ouder, en gezinnen worden kleiner. Ontwikkeling leidt mogelijk tot een verschuiving in de vraag - van eensgezinswoningen naar kleinere woningen.

Dit verband is niet duidelijk terug te zien in de data. Het aantal 65+ ers is wel flink gegroeid - maar de verhoudingen tussen de groepen blijft redelijk stabiel.

df <- 
  bevolking %>%
  select (date, kl_20_jaar, v20_40_jaar, v40_65_jaar, v65_80_jaar, gd_80_jaar) %>%
  melt (id = c('date')) %>%
  mutate (year = factor(year(date))) %>%
  mutate (variable = factor (variable) ) 

levels(df$variable) = c('< 20 jaar', 
                        '20-40 jaar', 
                        '40-65 jaar', 
                        '65-80 jaar', 
                        '> 80 jaar')
df %>% 
  mutate (value = value / 1000) %>%
  filter (date == as.Date("2020-01-01") |
            date == as.Date("2010-01-01") |
            date == as.Date("2000-01-01") |
            date == as.Date("1990-01-01") |
            date == as.Date("1980-01-01") ) %>%
  ggplot( aes( x = variable, y = value, fill = year) ) +
  geom_bar(color = "black",  stat="identity", position="dodge") +
  scale_fill_brewer()  + 
  coord_flip() +
  labs(
    title = "Bevolkingssamenstelling",
    y = "(x 1000)",
    x = ''
    ) + 
    guides(fill = guide_legend(reverse = TRUE))

bevolking %>% 
  select (date, kl_20_jaar, v20_40_jaar, v40_65_jaar, v65_80_jaar, gd_80_jaar) %>% 
  filter (date >= as.Date("1970-01-01")) %>% 
  melt (id = c('date')) %>% 
  mutate ( variable = factor(variable), 
           variable = factor(variable, levels = rev(levels(variable)))) %>% 
  mutate (value = value / 1000) %>% 
  ggplot ( aes( x = date, y = value, fill = variable) ) +
  geom_area(color = "darkgray", stat = "identity") +
  scale_fill_brewer(palette = "Set2", 
                    name = "Leeftijdsgroep",
                    direction = -1, labels = c("> 80 jaar", 
                                                          "65 > 80 jaar", 
                                                          "40 > 65 jaar",
                                                          "20 > 40 jaar",
                                                          "< 20 jaar"))  + 
  labs(
    title = "Bevolkingsgroei",
    subtitle = "1950 tot 2020",
    y = "Aantal (x 1000)",
    x = '',
    caption = "Bron: CBS, Bevolking; kerncijfers"
    )

De groei van de Nederlandse bevolking wordt in de laatste jaren steeds meer bepaald door migratie. Het geboorteoverschot (aantal nieuw geborenen -/- aantal sterfgevallen is nog nooit zo laag geweest). Mogelijk zien we hierdoor regionaal een sterkere vraag naar woningen in de steden.

bevolking %>%
  select (date, geboorteoverschot, migratiesaldo) %>% 
  filter (date >= as.Date("1970-01-01")) %>% 
  melt (id = c('date')) %>% 
  mutate (value = value / 1000) %>% 
  drop_na %>% 
  ggplot(aes (x = date, y=value, fill = variable)) +
  geom_bar(stat="identity") +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    title = "Bevolkingsgroei",
    y = "Aantal (x 1000)",
    x = ''
    )

Het CBS maakt regelmatig een prognose van de verwachte groei. Hier zien we dat deze trand zich de aankomende jaren voortzet. De groei van met name eenspersoonshuishoudens gaat de aankomende jaren nog zo door. De laatste prognoses zijn overigens van 2019.

aantal_hh <- 
  bevolking %>% 
  select (date,
            eenpersoonshuishoudens,
            meerpersoonshuishoudens) %>% 
  melt(id = c('date')) %>% 
  set_names ('date','type_huishouden', 'aantal')
  
prognose_hh <-
  prognose_huishoudens_2019 %>% 
  select (date,
            eenpersoonshuishoudens,
            meerpersoonshuishoudens) %>% 
  melt(id = c('date')) %>% 
  set_names ('date','type_huishouden', 'prognose') %>% 
  mutate(prognose = prognose / 1000)

aantal_hh %>% 
  full_join(prognose_hh, by = c("date", "type_huishouden")) %>% 
  melt (id = c('date','type_huishouden')) %>% 
  filter( date >= as.Date("1970-01-01")) %>%
  filter( date <= as.Date("2040-01-01")) %>% 
  drop_na() %>% 
  ggplot(aes (x = date)) +
  geom_line(aes(y=value,  color = type_huishouden, linetype = variable )) +
  scale_color_brewer(palette = "Accent") +
  labs(
    title = "Prognose huishoudens",
    y = "Aantal (x 1000)",
    x = '',
    colour = "Type huishoudens",
    linetype = ""
  )

5.1.2 Aanbod van woningen

Wanneer we een wat langere termijn trend kijken naar het aantal gerealiseerde nieuwbouwhuizen zien we dat het aantal nieuwbouwhuizen eigenlijk nu niet bijzonder hoog is. Ondanks dat er meer dan 70.000 huizen bijgebouwd wordt per jaar is dat vergeleken met de totale woningvoorraad eigenlijk niet bijzonder veel. Wel opvallend is dat de laatste jaren er duidelijk meer huizen onttrokken worden aan de woningvoorraad. Mogelijk geeft dit aan dat de ruimte voor nieuwbouw schaarser wordt. We zagen al dat de woningvoorraad lange termijn tred houdt met het aantal huishoudens.

voorraad_woningen %>%
  select(date,
         nieuwbouw,
         overige_toevoegingen,
         sloop,
         Overige_onttrekking,
         correctie
         ) %>%
  mutate (sloop = - sloop ) %>% 
  mutate (Overige_onttrekking = - Overige_onttrekking) %>% 
  mutate (correctie = - correctie) %>% 
  filter( date >= as.Date("1970-01-01")) %>% 
  melt (id = c('date')) %>% 
  mutate (value = value / 1000) %>% 
  ggplot(aes(x=date, y = value, fill = variable )) +
  geom_bar(stat="identity") +
  geom_hline(yintercept  = 0, color = "steelblue") +
  labs ( title = "Ontwikkelingen woningvoorraad",
         y = "Aantal (x 1000)",
         x = "Jaar") + 
    guides(fill = guide_legend(title = ''))

5.1.3 Samenhang Vraag en aanbod

#82235NED CBS Aantal Woningen
#Aantal huishoudens: NLHSTOTP Thomson Reuters

# Tidy data 
df <- 
  aantal_huishoudens_per_jaar %>%
  full_join(voorraad_woningen, by = c("date" = "date")) %>%
  select (date, 
          aantal_huishoudens, 
          voorraad) %>%
  set_names('date',
            'huishoudens',
            'woningvoorraad' ) %>%
  mutate(te_kort = (woningvoorraad - huishoudens)) %>%
  filter(date > as.Date("2000-01-01")) %>%
  melt(id = c('date')) 

# Line plot
lp<- 
  df %>%
  filter (variable != 'te_kort') %>%
  ggplot(aes (x = date)) +
  geom_line(aes(y=value, color = variable)) +
  labs(
    title = "Aantal huishoudens vs voorraad woningen",
    y = "Aantal (x 1000)",
    x = '',
    color = ''
    ) 

# Bar plot
bp <- 
  df %>%
  filter (variable == 'te_kort') %>%
  ggplot(aes (x = date)) +
  geom_bar(stat = "identity" , 
           aes(y = value, fill = factor(variable, labels = c("Te Kort woningen")))) +
  labs ( y = "Te kort (x 1000)",
         x = "",
         fill = "") +
  theme( axis.text.x = element_blank(),
         axis.ticks = element_blank()) 
  

ggarrange(lp, bp,  
          heights = c(2, 1),
          ncol = 1, nrow = 2, align = "v"
          )

In bovenstaande grafiek zien we dat het aantal woningen al decenia lang gelijke tred houdt met het aantal woningen. Er is daarbij een klein tekort zichtbaar. Dit is veel kleiner dan het gerapporteerde tekort van 331.000 woningen [CBS2020b]. De prognose van het aantal huishoudens geeft een voortzetting van deze trend aan.

Er is in de media veel aandacht voor of de bouwplannen wel realiseerbaar zijn. Mogelijk kan door regelgeving (PFAS, ruimtelijke ordening) de groei van het aantal huizen in de toekomst de groei van de bevolking niet bijhouden waardoor het tekort nog verder oploopt.

Dit zou verder onderzocht kunnen worden. In deze exploratieve studie wordt hier niet verder op ingegaan. De voorlopige conclusie is dan ook dat er een constante relatie is tussen vraag en aanbod op de huizenmarkt voor Nederland als geheel.

De voorzichtige conclusie is daarom dat er geen invloed is aangetoond van een verschuiving in vraag en aanbod op de prijsvorming in de huizenmarkt in Nederland.

5.1.4 Indicatoren Huizenprijzen

  • Gem. huizenprijzen (prijsindex)
  • Besteedbaar Inkomen
  • Particulier vermogen
  • Werkloosheid
  • Rentestanden
  • Financiering hypotheek of eigen geld
  • Fiscale behandeling, subsidies
  • Toegang tot krediet
  • Voorkeuren van kopers
  • Percentage huur / koop woningen
  • Consumentenvertrouwen

Vanwege de (structurele) te korten op de woningmarkt lijken huizen prijzen in belangrijke mate beinvloed te worden door de financieringsmogelijkheden van de kopers. De rente is historisch laag, de inkomens zijn gestegen en het eigen vermogen is gegroeit. De fiscale behandeling en regelgeving speelt hier een belangrijke rol, zowel ruimte scheppend - verruimen van de leenmogelijkheden van tweeverdieners, als beperkend - beperking aftrekbaarheid tot 100% financiering.

# Huizenprijzen

lp_prijzen <-
  prijsindex_woningen_lt %>% 
  filter(date >= as.Date('1995-01-01')) %>% 
  ggplot(aes(x = date)) +
  geom_line(aes( y = huisprijsindex), color = 'steelblue') +
  labs(
    title = "Huizenprijzen",
    subtitle = "Index (2015 = 100)",
    y = "",
    x = ""
  ) 


# GDP
lp_gdp <- 
  gdp_lt %>% 
  filter(date >= as.Date('1995-01-01')) %>% 
  ggplot( aes(x=date, y=gdp_index)) +
    geom_line(color = 'steelblue') +
    labs(
      title = "Ontwikkeling trend BBP (Bruto Binnenlands Product)",
      y = "",
      x = ''
      ) 

# Rente
lp_rente <-
  hypotheek_rente %>%
  filter(period >= as.Date("1995-01-01")) %>% 
  ggplot( aes(x=period, y=hypotheek_rente)) +
  geom_line(color = 'steelblue')  +
  labs(
    title = "Hypotheekrente",
    y = "gemiddelde rente",
    x = ''
    ) 


# Werkgelegenheid
lp_werkgelegenheid <-
  werkgelegenheid_lt %>% 
  filter (date >= as.Date('1995-01-01')) %>% 
  ggplot(aes(x = period, y = werkgelegenheid)) +
  geom_line(color = 'steelblue') +
  labs(
    title = "Werkgelegenheid",
    y = "(x 1000)",
    x = ''
  ) 

# Inkomen
lp_inkomen <-
  inkomen_lt %>% 
  filter (date >= as.Date('1995-01-01')) %>% 
  ggplot(aes(x = period, y = inkomen)) +
  geom_line(color = 'steelblue') +
  labs(
    title = "Inkomen per medewerker per kwartaal",
    y = "EUR",
    x = ''
  ) 

# Consumentenvertrouwen
lp_conf <-
  consumenten_vertrouwen_lt_q %>% 
  filter (date >= as.Date('1995-01-01')) %>% 
  mutate(pos = ifelse(consumer_conf >= 0,consumer_conf, 0)) %>% 
  mutate(neg = ifelse(consumer_conf < 0,consumer_conf, 0)) %>% 
  ggplot(aes(x = date)) +
  geom_area(aes(y = pos), fill = 'steelblue')  +
  geom_area(aes(y = neg), fill = 'indianred2')  +
  geom_hline (yintercept = 0, color = 'steelblue') +
  labs(
    title = "Consumenten vertrouwen indicator",
    y = "",
    x = ""
  ) 

ggarrange(lp_prijzen, lp_gdp, lp_rente,lp_inkomen, lp_werkgelegenheid,lp_conf,
          ncol = 2, nrow = 3, align = "v"
          )

Over de vermogens van huishoudens is er helaas niet een langere termijn tijdsreeks beschikbaar. Een veronderstelling van dit onderzoek was dat door het gewijzigdde beleid en de stijgende huizenprijzen het vermogen van huishoudens zou stijgen. Daardoor heeft men ook weer meer vermogen beschikbaar voor de aanschaf van woningen wat een prijsopdrijvend effect zou hebben. Dit is door gebrek en data helaas niet verder te verifieren.

df_vermogen <-
  vermogen %>% 
  filter(kenmerk_huishouden == "Particuliere huishoudens") %>%
  filter(vermogensbestanddeel == 'Vermogen exclusief eigen woning' )

inkomen %>%
  filter(kenmerk_huishouden == "Particuliere huishoudens")  %>% 
  select (date, besteedbaar_inkomen ) %>% 
  left_join( df_vermogen, by = (c('date')) ) %>% 
  select (date, besteedbaar_inkomen, mediaan_vermogen  ) %>% 
  melt (id = c('date')) %>%
  ggplot( aes(x=date, y=value, color = variable )) +
  geom_line() +
  theme_minimal() +
  labs(
    title = "Besteedbaar inkomen en vermogen huishoudens",
    y = "(x 1000)",
    x = ''
    ) +
  guides(color = guide_legend(title = '')) +
  scale_color_discrete(labels = c('Besteedbaar inkomen','Mediaan Vermogen'))

rm(df_vermogen)

5.1.5 Regressie Huizenprijzen

In deze fase voeren we een regressie analyse uit om te zien welke factoren de meest voorspellende waarde hebben voor de huizenprijzen.

Correlatie van huisprijzen

We vinden een zeer sterke correlatie tussen inkomen, werkgelegenheid, BBP en huizenprijzen. De verwachting is dat dit met name wordt veroorzaakt door de tijdswaarde. Over de jaren heen, stijgt het inkomen, daalt de rente, en stijgen de huizenprijzen.

# plot correlations
corMatrix <- round(cor(df_rm_huizenprijzen %>% 
                         select(-date)),5)

corCols <- c("huisprijsindex", 
             "werkgelegenheid",
             "BBP",
             "Cons.Vertrouwen", 
             "Inkomen")

colnames(corMatrix) <- corCols
rownames(corMatrix) <- corCols

corrplot(corMatrix,
            title = "Correlatie huizenprijzen",
            type = "lower", 
            # tl.cex = 0.8,
            tl.col = "black", 
            tl.srt = 45)

Lineair regressie model voor inkomen en huizenprijzen.

df_rm_huizenprijzen %>% 
  ggplot(aes(x = inkomenindex, y = huisprijsindex)) +
  geom_point() +
   geom_smooth(method = "lm")

# Inkomen vs huizenprijzen

lm_inkomen = lm(huisprijsindex~inkomenindex, 
                data = df_rm_huizenprijzen) 

tab_model(lm_inkomen)
  huisprijsindex
Predictors Estimates CI p
(Intercept) -49.04 -55.95 – -42.12 <0.001
inkomenindex 1.68 1.59 – 1.76 <0.001
Observations 234
R2 / R2 adjusted 0.868 / 0.867

Lineair regressie model voor bruto binnenlands product en huizenprijzen.

# GDP geeft de beste fit

df_rm_huizenprijzen %>% 
  drop_na() %>% 
  ggplot(aes(x = date, y = gdp_index)) +
  geom_line(color = 'steelblue') +
  labs(
    title = "gdp incremental",
    y = "",
    x = ''
  ) 

lm_gdp = lm(huisprijsindex~gdp_index, data = df_rm_huizenprijzen)

tab_model(lm_gdp)
  huisprijsindex
Predictors Estimates CI p
(Intercept) -122.09 -129.99 – -114.19 <0.001
gdp_index 0.38 0.37 – 0.40 <0.001
Observations 234
R2 / R2 adjusted 0.921 / 0.921 
df_rm_huizenprijzen %>% 
  ggplot(aes(x = gdp_index, y = huisprijsindex)) +
  geom_point() +
  geom_smooth(method = "lm")

df_rm_lm <- df_rm_huizenprijzen

# Residuals tonen wel nog een patroon - dus er zit meer in de data
df_rm_lm$predicted <- predict(lm_gdp)   # Save the predicted values
df_rm_lm$residuals <- residuals(lm_gdp) # Save the residual values

ggplot(df_rm_lm, aes(x = gdp_index, y = residuals )) +
  geom_point(color = 'red', alpha = 0.5) +
  labs(
    title = "Residuals BBP Lineair Model",
    y = "Residuals",
    x = 'BBP Index'
  )

De residuals na deze correlatie laten wel zien dat er nog andere factoren zijn die van belang zijn voor het voorspellen van de huizenprijzen.

Wanneer het model wordt uitgebreid met consumer confidence, kan marginaal nog een iets hogere correlatie worden gevonden. Toevoegen van werkgelegenheid of inkomen geeft geen betere resultaten.

lm_gdp2 = lm(huisprijsindex~gdp_index + consumer_conf, 
             data = df_rm_huizenprijzen) 

#Create a linear regression with a quadratic coefficient
tab_model(lm_gdp2) #Review the results
  huisprijsindex
Predictors Estimates CI p
(Intercept) -117.07 -124.83 – -109.31 <0.001
gdp_index 0.37 0.36 – 0.39 <0.001
consumer_conf -0.15 -0.21 – -0.09 <0.001
Observations 234
R2 / R2 adjusted 0.929 / 0.929 
#Plot the Cooks Distances.
ggplot(lm_gdp2, aes(seq_along(.cooksd), .cooksd)) +
  geom_col(color = 'steelblue') 

# Multiple regressie, wordt de relatie niet heel veel beter van
lm_gdp_wg = lm(huisprijsindex~gdp_index+werkgelegenheid, 
               data = df_rm_huizenprijzen) 
tab_model(lm_gdp_wg)
  huisprijsindex
Predictors Estimates CI p
(Intercept) -119.50 -158.10 – -80.91 <0.001
gdp_index 0.39 0.23 – 0.56 <0.001
werkgelegenheid -0.00 -0.02 – 0.02 0.893
Observations 234
R2 / R2 adjusted 0.921 / 0.921 
# Multiple regressie, wordt de relatie niet heel veel beter van
lm_gdp_ink = lm(huisprijsindex~gdp_index+inkomenindex, 
                data = df_rm_huizenprijzen) 
tab_model(lm_gdp_ink)
  huisprijsindex
Predictors Estimates CI p
(Intercept) -115.43 -127.04 – -103.82 <0.001
gdp_index 0.34 0.29 – 0.40 <0.001
inkomenindex 0.19 -0.05 – 0.43 0.125
Observations 234
R2 / R2 adjusted 0.922 / 0.921

5.1.6 Huizenprijzen voorspellen

Tot slot is voor een beter passend model nog gekeken naar MARS (multivariate adaptive regression splines [Friedman (1991a, 1991b)]. Deze methode splits de regressie in meerdere ‘splines’. Waardoor een betere fit gevonden kan worden door rekening te houden met een non-lineaire relatie over de tijd.

Het MARS model kiest zelf de meest relevante parameters. Het MARS model komt hierbij uit op Inkomen en GBP. Wanneer met dit model geprobeerd wordt om de prijsontwikkeling van de afgelopen 5 jaar te voorspellen geeft dit model geen betere match.

# Simple MARS model
df_train <- df_rm_huizenprijzen %>% 
  filter(date <= (max(df_rm_huizenprijzen$date) - years(5)))
df_test <- df_rm_huizenprijzen %>% 
  filter(date > (max(df_rm_huizenprijzen$date) - years(5)))


max(df_test$date)
## [1] "2020-10-01"
# Tuning the model
# Optimal number of defrees and number of variables to prune is calculated

# create a tuning grid 
hyper_grid <- expand.grid(
  degree = 1:3, 
  nprune = seq(2, 7)
  )


# for reproducibiity
set.seed(123)

# cross validated model
tuned_mars <- train(
  x = subset(df_train, select = -huisprijsindex),
  y = df_train$huisprijsindex,
  method = "earth",
  metric = "RMSE",
  trControl = trainControl(method = "cv", number = 10),
  tuneGrid = hyper_grid
)

# plot results
ggplot(tuned_mars)

# fit model
fit <- earth(huisprijsindex~., df_train, 
             degree = tuned_mars$bestTune$degree, nprune = tuned_mars$bestTune$nprune)
# summarize the fit
summary(fit)
## Call: earth(formula=huisprijsindex~., data=df_train,
##             degree=tuned_mars$bestTune$degree,
##             nprune=tuned_mars$bestTune$nprune)
## 
##                                                coefficients
## (Intercept)                                       74.515567
## h(434.381-gdp_index)                              -0.093637
## h(gdp_index-434.381)                               0.271870
## h(90.1389-inkomenindex)                           -1.365846
## h(inkomenindex-90.1389)                           -2.447051
## h(528.546-gdp_index) * h(90.1389-inkomenindex)     0.002257
## h(gdp_index-528.546) * h(90.1389-inkomenindex)     0.012713
## 
## Selected 7 of 7 terms, and 2 of 5 predictors (nprune=7)
## Termination condition: RSq changed by less than 0.001 at 7 terms
## Importance: inkomenindex, gdp_index, date-unused, werkgelegenheid-unused, ...
## Number of terms at each degree of interaction: 1 4 2
## GCV 2.212228    RSS 392.0749    GRSq 0.9978434    RSq 0.9981446
# summarize the importance of input variables
evimp(fit)
##              nsubsets   gcv    rss
## inkomenindex        6 100.0  100.0
## gdp_index           4   9.7    9.6
# make predictions
train_predictions <- predict(fit, df_train)
test_predictions <- predict(fit, df_test)

df_result <-
  tibble (
    date = df_train$date, 
    huisprijsindex = df_train$huisprijsindex,  
    predicted = train_predictions,
    variable = 'train'
  ) %>% 
  union_all( 
    tibble(
      date = df_test$date, 
      huisprijsindex = df_test$huisprijsindex,  
      predicted = test_predictions,
      variable = 'test')
)

df_result %>% 
  ggplot(aes(date, huisprijsindex)) +
  geom_point(size = 1) +
  geom_line(aes(y = predicted, color = variable), size = 1, alpha = 1) +
  labs(
    title = "Voorspelling huizenprijzen MARS Model",
    y = "huisprijsindex",
    x = '',
    color = 'dataset'
  )

5.1.7 Regionale verschillen

Er bestaan grote regionale verschillen in huisprijzen, en in de prijsonwikkeling over de tijd.

# Distribution
verkoopprijs_regio %>% 
  filter(date >= as.Date('2015-01-01')
            ) %>% 
  drop_na(verkoopprijs) %>% 
  mutate(verkoopprijs = verkoopprijs / 1000) %>% 
  ggplot(aes(x = verkoopprijs, color = '', fill = '') ) +
  geom_density(alpha = 0.2) +
  scale_fill_manual( values = c("steelblue")) +
  scale_color_manual( values = c("steelblue")) +
  facet_wrap(~jaar) +
  labs(title = "Distributie verkoopprijzen (2015-2020)",
        x = "Verkoopprijs (x 1000 eur)") + 
  theme(legend.position = "none")

# Create a thematic map
verkoopprijs_regio %>%
  filter(date >= as.Date('2015-01-01')
           ) %>% 
  mutate(verkoopprijs = verkoopprijs / 1000) %>% 
  ggplot() +
  geom_sf(aes(fill = verkoopprijs)) +
  facet_wrap(~jaar) +
  scale_fill_viridis_c(option = "inferno") +
  labs(title = "Gemiddelde verkoopprijzen woningen", 
       subtitle = "",
       fill = "(x 1000 eur)") +
  theme_void()

# Add house prices per region (Boxplots)

5.2 Aantal transacties

Het aantal transacties heeft inde crisis jaren 2009 - 2012 een duidelijke terugval laten zien. Sindsdien zien we een snele inhaalbeweging, en nog steeds zien we aanzienlijk hogere transactie volumes dan voorheen.

verkochte_woningen_per_type_per_kwartaal %>% 
  mutate(subtype = factor(subtype, levels=c('Vrijstaande woning',
                                  '2-onder-1-kapwoning',
                                  'Hoekwoning',
                                  'Tussenwoning',
                                  'Appartement',
                                  'Onbekend'))) %>% 
  ggplot(aes(x=date, y=verkochte_bestaande_woningen )) +
  geom_line(aes(color = subtype)) +
  facet_wrap(~ subtype, ncol = 2, scales = "free") +
  labs(
    title = "Verkopen Bestaande Koopwoningen per kwartaal",
    y = "",
    x = '',
    color = 'subtype'
    ) +
  theme(legend.position = "none")

Wanneer we kijken naar het aantal transacties valt verder de piek op in apartementen op in Q4 2020. Dit is te verklaren door de wijziging in de regelgeving voor beleggers. Particuliere beleggens op de woningmarkt moeten vanaf 2021 8% overdrachtsbelasting betalen. Dit heeft gezorgt voor een piek aan aankopen in het laatste kwartaal van 2020.

Voorspellen van aantal transacties

Ook hier wordt eerst gekeken naar de correlaties

corMatrix <- round(cor(df_rm_transactions %>% 
                         select(-date)),5)

corCols <- c("huisprijsindex", 
             "Inkomen",
             "werkgelegenheid",
             "BBP",
             "Cons.Vertrouwen", 
             "Aantal huishoudens",
             "Aantal transacties")

colnames(corMatrix) <- corCols
rownames(corMatrix) <- corCols

corrplot(corMatrix,
            title = "Correlatie transacties",
            type = "lower", 
            tl.cex = 0.8,
            tl.col = "black", 
            tl.srt = 45)

Seizoensinvloeden bij het aantal transacties

Om de seizoensinvloeden verder te analyseren heb ik een tijdreeksanalyse uitgevoerd op de data. De ACf (auto correlatie) plot toont dat meer recente data significant is, wat een trend aangeeft in de data. De PACF plot (partial auto-correlation function) toont dat recente data en seizoensinvloeden een goede indicatie geven voor het modeleren van de trend.

De uiteindelijke voorspelling laat duidelijk het seizoenspatroon zien. Verder wordt door dit model de recente trend doorgetrokken, waar eigenlijk de verwachting zou zijn dat op termijn de huizenmarkt terug keert naar een stabiel aantal transacties.

# 
# vw <- verkochte_woningen_per_type_per_kwartaal %>% 
#   group_by (yearqtr) %>% 
#   summarise(totaal = sum(verkochte_bestaande_woningen) / 1000 )

# Time series object
vwts <- ts(verkochte_woningen_per_kwartaal["totaal"], 
              start = c(1995, 1), frequency = 4)

v1 <- 
  vwts %>%
  autoplot() +
  geom_smooth() +
  labs(
    title = "Verkopen Bestaande Koopwoningen per kwartaal",
    y = "(x 1000)",
    x = ''
    ) 

v2 <- 
  vwts %>% 
  ggsubseriesplot() +
  labs(
    title = "",
    y = "(x 1000)",
    x = ''
    ) 

# Plot it
cowplot::plot_grid(v1, v2, ncol=1, rel_heights = c(1, 1.5))

# Lag plot of non-seasonal effects
gglagplot(vwts, do.lines = FALSE) +
  labs(
    title = "Seizoensinvloeden Verkopen Bestaande Koopwoningen",
    y = "",
    x = '',
    color = 'Quarter'
    ) 

# ACF plot
# Visualizes how much the most recent value of the series is correlated with past values of the series
o1 <- ggAcf(vwts, lag.max = 16) +
  labs(
    title = "ACF Plot seizoensinvloeden Verkopen Bestaande Koopwoningen",
    y = "",
    x = ""
    ) 

# PACF Plot
# Visualizes whether certain lags are good for modeling or not; useful for data with a seasonal pattern
o2 <- ggPacf(vwts, lag.max = 16) +
  labs(
    title = "PACF Plot seizoensinvloeden Verkopen Bestaande Koopwoningen",
    y = "",
    x = ""
    ) 

# Plot both
cowplot::plot_grid(o1, o2, ncol = 1)

# Take first difference of the plot + seseonal effects

vwts %>% diff(lag=4) %>% diff() %>% ggtsdisplay(main = "Eerst orde verschil en seizoensinvloed")

# Eventualy I got this model - which seems to fit beter then auto arima
fit <-
  Arima(vwts, order=c(0,1,4), seasonal=c(0,1,1))
checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,4)(0,1,1)[4]
## Q* = 3.6063, df = 3, p-value = 0.3072
## 
## Model df: 5.   Total lags used: 8
fit %>% 
  forecast(h=12) %>% 
  autoplot() +
  labs(
    title = "Verwachting verkopen bestaande koopwoningen",
    y = "",
    x = ""
    )

Voorspellen aantal transacties

Ook is het mogelijk om net als bij de huizenprijzen, een MARS model toe te passen op het aantal transacties.

Het aantal transactie geeft een redelijke prognose. In ieder geval een betere aansluiting dan bij de huizenprijzen. Het MARS model heeft op basis van deze dataset en de gekozen variabelen de periode en het aantal huishoudens aangewezen als de belangrijkse voorspellende factoren.

df_train <- df_rm_transactions %>% 
  filter(date <= (max(df_rm_transactions$date) - years(5)))
df_test <- df_rm_transactions %>% 
  filter(date > (max(df_rm_transactions$date) - years(5)))


# Tuning the model
# Optimal number of defrees and number of variables to prune is calculated

# create a tuning grid 
hyper_grid <- expand.grid(
  degree = 1:3, 
  nprune = seq(2, 10)
  )


# for reproducibiity
set.seed(123)

# cross validated model
tuned_mars <- train(
  x = subset(df_train, select = -totaal),
  y = df_train$totaal,
  method = "earth",
  metric = "RMSE",
  trControl = trainControl(method = "cv", number = 27),
  tuneGrid = hyper_grid
)

# plot results
ggplot(tuned_mars)

# fit model
fit <- earth(totaal~., df_train, 
             degree = tuned_mars$bestTune$degree, nprune = tuned_mars$bestTune$nprune)
# summarize the fit
summary(fit)
## Call: earth(formula=totaal~., data=df_train, degree=tuned_mars$bestTune$degree,
##             nprune=tuned_mars$bestTune$nprune)
## 
##                            coefficients
## (Intercept)                   66.751759
## h(date-14791)                  0.031093
## h(16161-date)                 -0.023629
## h(113.803-huisprijsindex)     -1.628349
## h(inkomen-7685.59)            -0.052716
## h(inkomen-8244.91)             0.039445
## h(10048.9-inkomen)            -0.011964
## h(aantal_huishoudens-7313)     0.128697
## h(7569-aantal_huishoudens)     0.263887
## h(aantal_huishoudens-7569)    -0.300173
## 
## Selected 10 of 17 terms, and 4 of 7 predictors (nprune=10)
## Termination condition: Reached nk 21
## Importance: date, aantal_huishoudens, huisprijsindex, ...
## Number of terms at each degree of interaction: 1 9 (additive model)
## GCV 8.73284    RSS 1154.034    GRSq 0.9091409    RSq 0.9276717
# summarize the importance of input variables
evimp(fit)
##                        nsubsets   gcv    rss
## date                          7 100.0  100.0
## aantal_huishoudens            7 100.0  100.0
## huisprijsindex                6  94.7   94.9
## werkgelegenheid-unused        4  92.0>  95.5>
## inkomen                       3  39.4   39.5
# make predictions
train_predictions <- predict(fit, df_train)
test_predictions <- predict(fit, df_test)

df_result <-
  tibble (
    date = df_train$date, 
    totaal = df_train$totaal,  
    predicted = train_predictions,
    variable = 'train'
  ) %>% 
  union_all( 
    tibble(
      date = df_test$date, 
      totaal = df_test$totaal,  
      predicted = test_predictions,
      variable = 'test')
)

df_result %>% 
  ggplot(aes(date, totaal)) +
  geom_point(size = 1) +
  geom_line(aes(y = predicted, color = variable), size = 1, alpha = 1) +
  labs(
    title = "Voorspelling aantal transacties MARS Model",
    y = "Aantal transacties",
    x = '',
    color = 'dataset'
  )

6 Conclusies

Deze is ook leuk:

https://www.youtube.com/watch?v=5UvrPO4eY5A&feature=emb_logo

Referenties

CBS. (2016). Indicatoren woningmarkt op groen. Retrieved from https://www.cbs.nl/nl-nl/achtergrond/2016/16/indicatoren-woningmarkt-op-groen

CBS. (2020a). Huizenmarkt in beeld. Retrieved from https://www.cbs.nl/nl-nl/visualisaties/huizenmarkt-in-beeld

CBS. (2020b). Monitor koopwoningmarkt. Retrieved from https://www.cbs.nl/nl-nl/visualisaties/monitor-koopwoningmarkt

Foreman, J. W. (2014). Data Smart, Using Data Science to Transform Information into Insight. John Wiley & sons, Inc.

Groot, S., Vogt, B., Van der Wiel, K., & Van Dijk, M. (2018). Oververhitting op de Nederlandse huizenmarkt? Retrieved from https://www.cpb.nl/publicatie/oververhitting-op-de-nederlandse-huizenmarkt/CPB-Achtergronddocument-1jun2018-Oververhitting-op-de-nederlandse-huizenmarkt.pdf

Henry, L., & Wickham, H. (2020). Purrr: Functional programming tools. Retrieved from https://CRAN.R-project.org/package=purrr

Hyndman, R. J., & Athanasopoulos, G. (2018). Forecasting: Principles and Practice. Retrieved from https://otexts.com/fpp2/

Kenneth Gopal, Gerard van Leeuwen, David Omtzigt, T. K. en M. S.-F. (2020). Prognose woningmarktmodel Socrates. Retrieved from https://www.abfresearch.nl/publicaties/rapportage-socrates-2019/

Kenneth Gopal, Gerard van Leeuwen, David Omtzigt, T. K. en M., & Stuart-Fox. (2019). Socrates - Scenarioverkenningen van de woningmarkt in 2030. Retrieved from https://www.rijksoverheid.nl/documenten/rapporten/2019/05/29/socrates-2019---scenarioverkenningen-van-de-woningmarkt-in-2030

Léon Groenemeijer, Kenneth Gopal, D. O. &. G. van L. (2020). Vooruitzichten bevolking, huishoudens en woningmarkt, Prognose en Scenario’s 2020-2035. Retrieved from https://www.rijksoverheid.nl/documenten/rapporten/2020/06/12/vooruitzichten-bevolking-huishoudens-en-woningmarkt-prognose-en-scenarios-2020-2035

MacKay, R. J., & Oldford, R. W. (2000). Scientific method, statistical method and the speed of light. Statistical Science, 15(3), 254–278. https://doi.org/10.1214/ss/1009212817

Ministerie van Binnenlandse Zaken en Koninkrijksrelaties. (2021). Datawonen. Retrieved from https://datawonen.nl/dashboard/dashboard/koopprijs

Müller, K., & Wickham, H. (2021). Tibble: Simple data frames. Retrieved from https://CRAN.R-project.org/package=tibble

R Core Team. (2021). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from https://www.R-project.org/

Research, A. (2020). Primos - Bevolkingsprognose, Rapportage Primos. Retrieved from https://www.abfresearch.nl/producten/prognoses/primos-bevolkingsprognose/

Vogt, B., Kalara, N., & Voogt, B. (2018). How do the Dutch Finance their Own House? Retrieved from https://www.cpb.nl/sites/default/files/omnidownload/CPB-Background-Document-June2018-How-do-the-dutch-finance-their-own-house.pdf

Wickham, H. (2016). Ggplot2: Elegant graphics for data analysis. Springer-Verlag New York. Retrieved from https://ggplot2.tidyverse.org

Wickham, H. (2019a). Stringr: Simple, consistent wrappers for common string operations. Retrieved from https://CRAN.R-project.org/package=stringr

Wickham, H. (2019b). Tidyverse: Easily install and load the tidyverse. Retrieved from https://CRAN.R-project.org/package=tidyverse

Wickham, H. (2021a). Forcats: Tools for working with categorical variables (factors). Retrieved from https://CRAN.R-project.org/package=forcats

Wickham, H. (2021b). Tidyr: Tidy messy data. Retrieved from https://CRAN.R-project.org/package=tidyr

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., … Yutani, H. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686

Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., … Dunnington, D. (2020). Ggplot2: Create elegant data visualisations using the grammar of graphics. Retrieved from https://CRAN.R-project.org/package=ggplot2

Wickham, H., François, R., Henry, L., & Müller, K. (2021). Dplyr: A grammar of data manipulation. Retrieved from https://CRAN.R-project.org/package=dplyr

Wickham, H., & Hester, J. (2020). Readr: Read rectangular text data. Retrieved from https://CRAN.R-project.org/package=readr